#include <istream>
#include <vector>
#include <string>
#include <fstream>
#include <iostream>
#include <algorithm>
#include <map>
#include <pthread.h>
#include <boost/format.hpp>
#include <stdlib.h>

#include "bcpp_tree.h"
#include "portable_timer.h"
#include "perlutil.h"
#include "matrix.h"
#include "get_arg.h"
#define DEBUG_LEVEL 2
#include "debug_func.h"

using std::vector;
using std::ostringstream;

using std::string;
using std::fill;
using std::ofstream;
using std::cout;
using std::cin;
using std::cerr;
using std::ifstream;
using std::map;
using std::pair;
using std::istream;
using std::max;
using std::min;



//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	print and analyse orthologs

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
bool open_filename(ofstream& ofs, const string& file_name)
{
	ofs.open(file_name.c_str());
	if (ofs)
		return true;
	cerr << "Error\nCould not open " << file_name << "\n";
	return false;
}
bool print_gene_set(ofstream& ofs, const string& file_name, const string& header, const bcpp_tree::t_ortholog_data& data)
{
	if (!data.ensg_sets.size())
		return true;

	if (!ofs.is_open() &&
		!open_filename(ofs, file_name))
		return false;

//	cerr << boost::format("%|2$+3| sets from %|1$-10|:\n\t%|3$+5| tax1 genes\n\t%|4$+5| tax2 genes\n")
//	% file_name % data.ensg_sets.size() % data.cnt_tax1 % data.cnt_tax2;
	ofs << header;
	for (unsigned i = 0; i < data.ensg_sets.size(); ++i)
	{
		for (unsigned j = 0 ; j < data.ensg_sets[i].size(); ++j)
			ofs << data.ensg_sets[i][j].gene << ":"
				<< data.ensg_sets[i][j].len << "\t";
		ofs << "\n";
	}
	ofs << std::endl;
	return true;
}


struct t_outputstreams
{

	t_outputstreams(const string& prefix) : file_prefix(prefix)
	{
		open();

		ofs_all_nodes	<< "To leaf\tTo root\tGene id\tNode type\n"
						<< string(80, '-') << "\n";

		ofs_full_tree	<< "Branch length\tTo root\tTo leaf\tGene id\tNode type\n"
						<< string(80, '-') << "\n";
	}
	string	 file_prefix;
	ofstream timing;
	ofstream ortho_1_to_1;
	ofstream ortho_1_to_m;
	ofstream ortho_m_to_1;
	ofstream ortho_m_to_m;
	ofstream orphaned_1_to_1;
	ofstream orphaned_1_to_m;
	ofstream orphaned_m_to_1;
	ofstream orphaned_m_to_m;
	ofstream orphaned_tax1;
	ofstream orphaned_tax2;
	ofstream ofs_full_tree;
	ofstream newhampshire_tree;
	ofstream ofs_all_nodes;
	bool open();
	bool printout_analysis(	bcpp_tree::t_tree&		tree,
						const string&				header);
};


bool t_outputstreams::printout_analysis(bcpp_tree::t_tree& tree, const string& header)
{
//	cerr << "\n" << string(80, '*') << "\n\n";

	// dump all nodes
	{
		ofs_all_nodes << header << std::endl;
		bcpp_tree::do_dump_nodes dump_nodes(ofs_all_nodes);
		tree.for_each(dump_nodes);
		ofs_all_nodes << std::endl;
	}

	// print full tree
	{

		ofs_full_tree << header << std::endl;
		bcpp_tree::do_print_tree_full print_tree(ofs_full_tree);
		tree.for_each(print_tree);
		ofs_full_tree << std::endl;


		newhampshire_tree << header << std::endl;
		bcpp_tree::do_newhampshire_format_str nh_str(newhampshire_tree);
		tree.for_each(nh_str);
		newhampshire_tree << std::endl;
	}

	// print gene sets
	{
		bcpp_tree::do_collect_orthologs data;
		tree.for_each(data);
		if (!print_gene_set(ortho_1_to_1	, file_prefix + "ortho_1_to_1.data"     ,
													header, data.ortho_1_to_1		) ||
			!print_gene_set(ortho_1_to_m	, file_prefix + "ortho_1_to_m.data"     ,
													header, data.ortho_1_to_m   	) ||
			!print_gene_set(ortho_m_to_1	, file_prefix + "ortho_m_to_1.data"     ,
													header, data.ortho_m_to_1   	) ||
			!print_gene_set(ortho_m_to_m	, file_prefix + "ortho_m_to_m.data"     ,
													header, data.ortho_m_to_m   	) ||
			!print_gene_set(orphaned_1_to_1	, file_prefix + "orphaned_1_to_1.data"  ,
													header, data.orphaned_1_to_1	) ||
			!print_gene_set(orphaned_1_to_m	, file_prefix + "orphaned_1_to_m.data"  ,
													header, data.orphaned_1_to_m	) ||
			!print_gene_set(orphaned_m_to_1	, file_prefix + "orphaned_m_to_1.data"  ,
													header, data.orphaned_m_to_1	) ||
			!print_gene_set(orphaned_m_to_m	, file_prefix + "orphaned_m_to_m.data"  ,
													header, data.orphaned_m_to_m	) ||
			!print_gene_set(orphaned_tax1	, file_prefix + "orphaned_tax1.data"    ,
													header, data.orphaned_tax1  	) ||
			!print_gene_set(orphaned_tax2	, file_prefix + "orphaned_tax2.data"    ,
													header, data.orphaned_tax2  	))
		return false;
	}
	return true;

}



bool t_outputstreams::open()
{
	return(
			open_filename(ofs_full_tree, 		file_prefix + "full_topology.tree"	)	&&
			open_filename(newhampshire_tree,	file_prefix + "new_hampshire.tree"	)	&&
			open_filename(ofs_all_nodes, 		file_prefix + "all_nodes.list"		)	&&
			open_filename(timing, 				file_prefix + "timing"		)

		);
}


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
int main (int argc, char *argv[])
{
	// we use cout to print progress so flush per call
	cout.setf(std::ios_base::unitbuf);

	std::ios::sync_with_stdio(false);

	t_program_args args;
	if (!process_command_line_options(argc, argv, args))
		return 1;



	string tree_line;
	if (!std::getline(cin, tree_line))
	{
		cerr << "Error:\n\tCould not retrieve tree string\n";
		return 1;
	}

	bool success;
	string error;

	//
	//	read data from input file
	//
	t_matrix		matrix;
	vector<string> 	names;
	string			cluster_name;
	unsigned		cnt_missing_values;
	unsigned		cnt_duplicate_values;
	t_matrix_bool	missing_values;

	ifstream infile;
	{
		infile.open(args.ks_file.c_str());
		if (!infile)
		{
			std::cerr << "Error opening [" << args.ks_file << "] \n";
			return 1;
		}
	}

	unsigned		cnt_pairs =  read_matrix_from_stream_in_fasta_pairs_format
												(infile, matrix,
															missing_values,
															names,
															cluster_name,
															cnt_missing_values,
															cnt_duplicate_values);


	t_outputstreams strms (args.out_file_prefix);

	bcpp_tree::t_tree paralog_tree;
	boost::tie(success, error) =
			load_tree_from_newhampshire_format_str(tree_line, paralog_tree);

	// get sum squares for upgma
	bcpp_tree::do_sum_sqrs2 analyse_tree(matrix);
	paralog_tree.for_each(analyse_tree);
	double sum_sqrs, avg_percent_stddev;
	analyse_tree.results(sum_sqrs, avg_percent_stddev);

	strms.timing
			<< 	"sum_sqrs = "				<<  sum_sqrs
			<<	"avg_percent_stddev = "     <<  avg_percent_stddev << "\n";

	// set type for each node
	bcpp_tree::do_set_node_types set_node_types(args.taxon1, args.taxon2);
	paralog_tree.for_each(set_node_types);

	if (!strms.printout_analysis(paralog_tree, "tree"))
		return 0;


	cerr << "\nAll finished!!\n\n";

	return 0;
}

